       Program nozzle
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C    A Program of quasi-one dimensional nozzle flows
C    Written by Dong Gang, May, 2003
C
C    Based on section 7.3&7.4, chapter 7, Computational Fluid 
C    Dynamics, The Basics with Applications,
C    John D. Anderson, JR. McGraw-Hill, 2002, 4
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      parameter (n=31)
      dimension x(n),a(n),r(n),t(n),u(n),p(n),
     +          am(n),fm(n),rm(n),um(n),tm(n)
      dimension delt(n),drdt(n),dudt(n),dtdt(n),drdtm(n),dudtm(n),
     +          dtdtm(n),drdtav(n),dudtav(n),dtdtav(n)
      dimension rt(10000),tt(10000),pt(10000),amt(10000)
      gamma=1.4
      open(11,file='output\output.dat',status='unknown')
      open(12,file='output\outputa.dat',status='unknown')
      write(*,*)'please input the nstop'
      read (*,*)nstop
5     write(*,*)'Subsonic-Supersonic/Purely subsonic(0/1)?'
      read (*,*)icase
C
C    Initial conditions specified
C
      if (icase.eq.0) then 
        do i=1,n
          x(i)=(i-1)/((n-1)/3.)
          a(i)=1.+2.2*(x(i)-1.5)*(x(i)-1.5)
          r(i)=1.-0.3146*x(i)
          t(i)=1.-0.2314*x(i)
          u(i)=(0.1+1.09*x(i))*(t(i)**0.5)
        end do
      elseif (icase.eq.1) then
        write(*,*)'In purely subsonic case, input the exit pressure,'
        write(*,*)'the range is from 0.528P0 to P0'
        read(*,*) pe
        do i=1,n
          x(i)=(i-1)/((n-1)/3.)
          a(i)=1.+0.2223*(x(i)-1.5)*(x(i)-1.5)
          if (x(i).le.1.5) a(i)=1.+2.2*(x(i)-1.5)*(x(i)-1.5)         
          r(i)=1.-0.023*x(i)
          t(i)=1.-0.009333*x(i)
          u(i)=0.05+0.11*x(i)
        end do
      else
        write(*,*)'There is no such case, repeat now!'
        goto 5
      end if          
C
C    Time stepping and looping
C
      deltx=3./real(n-1)
      cfl=0.5
      do ii=1,nstop
      dt=1.
      do i=1,n
        as=sqrt(t(i))
        delt(i)=cfl*deltx/(as+u(i))
        if (dt.gt.delt(i)) dt=delt(i)
      end do
C
C    Get solutions by MacCormack's explicit technique
C
      do i=2,n-1
      drdt(i)=(-r(i)*(u(i+1)-u(i))-r(i)*u(i)*(log(a(i+1))-log(a(i)))
     +        -u(i)*(r(i+1)-r(i)))/deltx
      dudt(i)=(-u(i)*(u(i+1)-u(i))-1./gamma*(t(i+1)-t(i)+
     +        t(i)/r(i)*(r(i+1)-r(i))))/deltx
      dtdt(i)=(-u(i)*(t(i+1)-t(i))-(gamma-1.)*t(i)*(u(i+1)-u(i)+
     +        u(i)*(log(a(i+1))-log(a(i)))))/deltx
C
      rm(1)=r(1)
      um(1)=u(1)
      tm(1)=t(1)
      rm(i)=r(i)+drdt(i)*dt
      um(i)=u(i)+dudt(i)*dt
      tm(i)=t(i)+dtdt(i)*dt
C
      drdtm(i)=(-rm(i)*(um(i)-um(i-1))-rm(i)*um(i)*(log(a(i))-
     +         log(a(i-1)))-um(i)*(rm(i)-rm(i-1)))/deltx
      dudtm(i)=(-um(i)*(um(i)-um(i-1))-1./gamma*(tm(i)-tm(i-1)
     +         +tm(i)/rm(i)*(rm(i)-rm(i-1))))/deltx
      dtdtm(i)=(-um(i)*(tm(i)-tm(i-1))-(gamma-1.)*tm(i)*
     +         (um(i)-um(i-1)+um(i)*(log(a(i))-log(a(i-1)))))/deltx
C
      drdtav(i)=(drdt(i)+drdtm(i))/2.
      dudtav(i)=(dudt(i)+dudtm(i))/2.
      dtdtav(i)=(dtdt(i)+dtdtm(i))/2.
C
      r(i)=r(i)+drdtav(i)*dt
      u(i)=u(i)+dudtav(i)*dt
      t(i)=t(i)+dtdtav(i)*dt
      p(i)=r(i)*t(i)
      am(i)=u(i)/sqrt(t(i))
      fm(i)=r(i)*a(i)*u(i)
        if (i.eq.15) then
          rt(ii)=r(i)
          tt(ii)=t(i)
          pt(ii)=p(i)
          amt(ii)=am(i)
        end if
      end do
C
C    Boundary conditions at inflow (n=1)
C
      u(1)=2.*u(2)-u(3)
      r(1)=1.
      t(1)=1.
      p(1)=r(1)*t(1)
      am(1)=u(1)/sqrt(t(1))
      fm(1)=r(1)*a(1)*u(1)
C
C    Boundary conditions at outflow (n=n)
C
      if (icase.eq.0) then
        u(n)=2.*u(n-1)-u(n-2)
        r(n)=2.*r(n-1)-r(n-2)
        t(n)=2.*t(n-1)-t(n-2)
        p(n)=r(n)*t(n)
        am(n)=u(n)/sqrt(t(n))
        fm(n)=r(n)*a(n)*u(n)
      else
        p(n)=pe
        u(n)=2.*u(n-1)-u(n-2)
        t(n)=2.*t(n-1)-t(n-2)
        r(n)=p(n)/t(n)
        am(n)=u(n)/sqrt(t(n))
        fm(n)=r(n)*a(n)*u(n)
      end if       
C
      end do
C
C    Print the results
C
      do i=1,n
        write(11,100)x(i),a(i),r(i),u(i),t(i),p(i),am(i),fm(i)
      end do
      do ii=1,nstop
        write(12,101)ii,rt(ii),tt(ii),pt(ii),amt(ii)
      end do
100   format(f5.2,7f8.3)
101   format(i3,4f8.3)
      end